#include "AnaDef.h"
#include "XGLUtils.h"
#include "style.h"

const Double_t EPSILON = 1e-12;


void normHist(TH1* hh)
{
  for(Int_t ii=1; ii<=hh->GetNbinsX(); ii++){
    const Double_t cc = hh->GetBinContent(ii);
    const Double_t err = hh->GetBinError(ii);
    const Double_t bw = hh->GetBinWidth(ii);

    hh->SetBinContent(ii, cc/bw);
    hh->SetBinError(ii, err/bw);
  }
}


Int_t IDx(const Int_t grid)
{
  const Int_t idx = grid/10;
  if(idx!=1 && idx!=2){
    printf("bad idx %d %d\n", grid, idx); exit(1);
  }
  return idx;
}

Int_t IDy(const Int_t grid)
{
  const Int_t idy= grid%10;
  if(idy!=6 && idy!=7){
    printf("bad idy %d %d\n", grid, idy); exit(1);
  }
  return idy;
}

int main(int argc, char * argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=1){
    printf("argc!=1\n");return 1;
  }

  TH1::AddDirectory(kFALSE);

  const TString tuneName[]={"Perugia0", "Perugia2011", "Perugia2010NoCR"};
  const TString pname[]={"p","pi","K"};
  const TString jpt[]={"5-10","10-15","15-20","incl"};
  const Int_t gids[]={16,17,26,27};

  TFile *fin = 0x0;

  for(Int_t itune=0; itune<(Int_t)(sizeof(tuneName)/sizeof(TString)); itune++){
    for(Int_t ijp=0; ijp<(Int_t)(sizeof(jpt)/sizeof(TString)); ijp++){

      printf("\n\n itune %d ijp %d\n", itune, ijp);
      if(jpt[ijp]!="incl"){
        fin = new TFile("theory.root");
      }
      else{
        fin = new TFile(Form("outMergePythiaFastJet_%s_R04.root",tuneName[itune].Data()));
      }
      if(!fin->IsOpen()){
        printf("fin not open %s\n", fin->GetName()); exit(1);
      }
      fin->ls();
      printf("\n\n");


      for(Int_t id=0; id<(Int_t)(sizeof(gids)/sizeof(Int_t)); id++){
        if(jpt[ijp]=="incl" && gids[id]!=17){
          continue;
        }

        TGraphAsymmErrors * grstat[]={0x0,0x0,0x0,0x0};

        TString stmp;
        //no electron
        for(Int_t itype=0; itype<3; itype++){
          TH1F * hin = 0x0;
          if(jpt[ijp]!="incl"){
            stmp=Form("%sGen_xg_%s", tuneName[itune].Data(), pname[itype].Data());
            TDirectoryFile * dd = (TDirectoryFile*) fin->Get(stmp);
            if(!dd){
              printf("no directory %s\n", stmp.Data()); exit(1);
            }
            dd->ls();

            TString tmpjp = jpt[ijp];
            tmpjp.ReplaceAll("5-10","05-10");
            tmpjp.ReplaceAll("-","_");
            
            TString xvar;
            if(IDx(gids[id])==1){
              xvar="TrackPt";
            }
            else if(IDx(gids[id])==2){
              xvar="Z";
            }
            else{
              printf("bad gids id %d %d\n", id, gids[id]); exit(1);
            }
            
            stmp=Form("fh1FF%sGen_xg_%s_%s", xvar.Data(), pname[itype].Data(), tmpjp.Data());
            hin = (TH1F*)dd->Get(stmp);

            if(!hin){
              printf("hist not found %s\n",stmp.Data()); exit(1);
            }

            //normalization only for pT hist
            if(IDx(gids[id])==1){
              style::ToNaturalScale(hin);
              normHist(hin);
            }

          }
          else{
            stmp = Form("hTrackPtGen_%s", pname[itype].Data());
            hin = (TH1F*) fin->Get(stmp);
            hin->RebinX(10);
            if(!hin){
              printf("hist not found %s\n",stmp.Data()); exit(1);
            }
          }

          printf("histgot %s\n", stmp.Data());

          //static TGraphAsymmErrors * HistToGraph(const TString grname, const TH1 *hh, const Double_t thres=0, const TH1 *herr=0x0, const Double_t xmin=-1e10, const Double_t xmax=1e10, const Bool_t klogx=kFALSE);
          //important to have klogx=1
          grstat[itype]=XGLUtils::HistToGraph(Form("grstat%d",itype), hin, 0, 0x0, -1e10, 1e10, IDx(gids[id])==1 && jpt[ijp]!="incl");
        }

        //to-pion ratio
        if(IDy(gids[id])==7){
          for(Int_t itype=0; itype<3; itype++){
            if(itype==1){
              continue;
            }

            for(Int_t ipt=0; ipt<grstat[itype]->GetN(); ipt++){
              const Double_t ypar = grstat[itype]->GetY()[ipt];
              if(ypar<EPSILON){
                continue;
              }
              const Double_t epar = grstat[itype]->GetEYlow()[ipt];

              const Double_t xtarget = grstat[itype]->GetX()[ipt];
              if( jpt[ijp]=="incl" && xtarget > 30){
                continue;
              }

              Int_t ipion = -999;
              for(Int_t itmp=0; itmp<grstat[1]->GetN(); itmp++){
                if(fabs( xtarget - grstat[1]->GetX()[itmp] )<EPSILON){
                  ipion = itmp;
                }
              }
              if(ipion<0){
                printf("ipion not found %e\n", xtarget); 
                for(Int_t itmp=0; itmp<grstat[1]->GetN(); itmp++){
                  printf("%d %e\n", itmp, grstat[1]->GetX()[itmp]);
                }
                exit(1);
              }
              //bad bad have to search
              const Double_t ypion = grstat[1]->GetY()[ipion];

              if(fabs(grstat[itype]->GetX()[ipt] - grstat[1]->GetX()[ipion])>EPSILON){
                printf("bad itype %d ipt %d %e pion %e\n", itype, ipt, grstat[itype]->GetX()[ipt], grstat[1]->GetX()[ipion]); exit(1);
              }
              const Double_t epion = grstat[1]->GetEYlow()[ipion];

              //====

              const Double_t newy = ypar/(ypion+EPSILON);
              const Double_t newE = newy*TMath::Sqrt( TMath::Power(epion/(ypion+EPSILON),2) + TMath::Power(epar/ypar,2) );//can NOT use sigma~sqrt(ypar), sqrt(ypion) because these are normalized yields!!

              grstat[itype]->SetPoint(ipt, grstat[itype]->GetX()[ipt], newy);
              grstat[itype]->SetPointError(ipt, grstat[itype]->GetEXlow()[ipt], grstat[itype]->GetEXlow()[ipt], newE, newE);
            }
          }
        }

        //electron dummy
        grstat[3]=new TGraphAsymmErrors; grstat[3]->SetName("grstat3");

        //grsys
        TGraphAsymmErrors * grsys[]={0x0,0x0,0x0,0x0};
        for(Int_t itype=0; itype<4; itype++){
          grsys[itype]=(TGraphAsymmErrors*)grstat[itype]->Clone(Form("grsys%d",itype));
          for(Int_t ipt=0; ipt<grsys[itype]->GetN(); ipt++){
            grsys[itype]->SetPointEYlow(ipt,0);
            grsys[itype]->SetPointEYhigh(ipt,0);
          }
        }

        TFile *fout=new TFile(Form("outTheoryCombined_%s_JPT%s_GrID%d.root", tuneName[itune].Data(), jpt[ijp].Data(), gids[id]),"recreate");
        for(Int_t itype=0; itype<4; itype++){
          grstat[itype]->Write();
          grsys[itype]->Write();
        }
        fout->Save();
        fout->Close();
        delete fout;
      }
    }
  }

  printf("hi convertTheorydonedone\n");
  return 0;
}
